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We report a Phase Switch Monte Carlo (PSMC) method study of the freezing line of the Lennard- 
Jones (LJ) fluid. Our work generalizes to soft potentials the original application of the method to 
hard sphere freezing, and builds on a previous PSMC study of the LJ system by Errington (J. Chem. 
Phys. 120, 3130 (2004)). The latter work is extended by tracing a large section of the Lennard- 
^ , Jones freezing curve, the results for which we compare to a previous Gibbs-Duhem integration study. 

O ' Additionally we provide new background regarding the statistical mechanical basis of the PSMC 

1/'^ ' method and extensive implementation details. 

■ PACS numbers; 64.70Fx, 68.35.Rh 
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I. INTRODUCTION AND BACKGROUND 

O 

^ The freezing of a disordered fluid into an ordered crystalline solid undoubtedly represents one of the most spectacular 

examples of thermodynamics in action. But despite its ubiquity and familiarity, key aspects of the phenomenon 
remain poorly understood across a variety of systems [1]. Principal among the reasons for this is the difficulty of 
obtaining accurate simulation estimates for fluid-solid phase coexistence properties. In this section we outline the most 
C/3 ] commonly used contemporary approaches for tracing freezing boundaries, identify their key strengths and weaknesses, 
and highlight recent developments in the search for improvements. 

The staple method for locating fluid-solid coexistence is thermodynamic integration (TI) [1~3]. Here, the free 
energy of each phase (fluid (F) and crystalline solid (CS)) is computed for states covering a range of densities, using 
integration methods which connect their thermodynamic properties with those of reference states of known free energy. 
fH The two branches of the free-energy are then matched to determine the freezing parameters. While there is much to 
Q ' commend TI (principally its simplicity) it turns out to be less than ideal in a number of respects [4] . These have been 
O ', documented in detail elsewhere (see eg. [5]); but briefly, the main difficulty is one of identifying a suitable reference 
state together with a reversible phase space pathway by which it can be reached; a poor choice of integration path 
' may encounter singularities -both real and artificial [6] - which compromise the calculation. Additionally, corrections 
J> . may be required to allow for the fact that the path does not quite reach the idealized reference state [2, 7]. Moreover, 
C ' the method focuses on the absolute free energies of the phases [8], rather than the quantity of interest -their free 
^ : energy difference. Finally, there currently exists no reliable and comprehensive method for estimating errors on phase 
] ' boundaries computed via TI. 

^ , Obtaining a whole F-CS phase boundary using TI is a computationally challenging task. However, knowledge 
ly-^ ■ of one point on the F-CS coexistence boundary (obtained, for instance, by using TI) can be used to bootstrap a 
Gibbs-Duhem integration (GDI) scheme to trace the entire curve without further calculation of free energies. The 
. GDI method [9-11] exploits the generalized Clausius-Clapeyron equation to express the slope of the phase boundary 
' entirely in terms of single-phase averages. This is clearly a virtue since it avoids the need for two-phase sampling. 
, However, without any 'reconnection' of the two configuration-spaces at subsequent simulation state points, the GDI 
I ' approach offers no feedback on integration errors. Since there will generally exist a band of metastable states on each 
' ^ ] side of the phase boundary, it is possible for the integration to wander significantly from the true boundary with no 
indication that anything is wrong. As a result, the calculation of meaningful uncertainties is problematic. GDI has 
been used in a number of studies, most notably in the context of freezing of hard and soft spheres [12]. 

More recently, attention has shifted to developing methods for tackling the problem of locating F-CS coexistence 
via direct measurements of free energy differences. To do so, one must construct a reversible sampling path between 
the phases. The obvious "physical" path, traversing the region where both phases coexist -whilst practicable when 
' the two phases share the same symmetry (see eg. ref. [13]) -turns out to be computationally problematic in the F-CS 
5^ , context. The main reasons for this are the large degree of metastability of the two phases, the extended timescale for 
the crystallization process, and the tendency of any crystal formed to exhibit vacancies, interstitials and other defects. 
Accordingly, recent effort has focussed on trying to identify alternative inter-phase routes for the F-CS problem. 

One such alternative route, called constrained fluid-A integration has been proposed by Grochla[14] and extended 
by Elke et al [15]. It involves the controlled transformation of the fluid phase to the solid phase in a series of stages. 
During the initial stage the fluid is gradually changed to a weakly interacting fluid by reducing the strength of 
interparticle interactions. At the next stage, a set of Gaussian potential wells of prescribed width are switched on at 
the sites of a crystalline lattice of the appropriate symmetry. Simultaneously, the volume of the system is gradually 
changed from a value typical of the fluid to one typical of the solid. The particles of the weakly interacting fluid are 
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then allowed to diffuse to the lattice sites, where they arc captured by the Gaussian wells. The final stage involves 
gradually restoring the particle interactions to full strength, whilst simultaneously switching off the Gaussian potential 
wells. Integration of free energy derivatives is used to estimate the free energy difference along the path. 

The constrained fluid-A integration method was tested by application to the Lennard-Jones fluid where it was used 
to determine two coexistence state points. This test led the author to conclude that the method "cannot be said 
to be computationally efficient" [14], at least in its present form. This finding presumably reflects the rather long 
and fragmented nature of the inter-phase path. Although each of the stages of the path was initially reported to 
be reversible [15, 16], further application of the method [16] found this to be contingent on the correct tuning (via 
trial and error) of four separate constants relating to the characteristics of the Gaussian potential wells. Indeed, it is 
not clear to us that one can generally expect a-priori that the stated path will be trouble-free; for example, it seems 
hard to rule out that the transformation of the original fluid to a weakly interacting fluid might intersect a first order 
(pseudo) liquid-vapor transition. Furthermore we speculate that for large systems the reliance on particle diffusion 
to translate particles to unoccupied lattice sites may become problematic. 

Very recently, Mastny and de Pablo [17] have proposed a method based on measurements of the density of states, 
which aims to directly estimate the free energy difference of F and CS phases. The rationale for the method is its 
authors' assertion that: "To connect the free energy of the solid and liquid phases, there must exist a portion of energy 
and volume space that can be simultaneously sampled by both solid and liriuid phases". To exploit this premise, a 
series of simulations of the Lennard-Jones system were performed, each restricted to -and sampling the configuration 
space of- one of a set of overlapping "windows" in the space of the total energy E and volume V. Successive windows 
were positioned appropriately such as to form a path linking the two well-separated regions of E, V associated with 
the typical configurations of the respective phases. Along this path, a central range of E, V was indeed reported to be 
found for which both liquid and solid phase configurations could be sampled. A simple average [18] of the solid phase 
and liquid phase density of states was accumulated in the central range and joined for continuity to the measured 
forms of the liquid and solid density of states on either side. 

In our opinion, the existence of a range of E and V that can be sampled by both liquid and solid phases is 
insufficient to connect the respective branches of the density of states because these are (like their underlying sets 
of characteristic configurations), fundamentally distinct. Instead, for a proper connection, the phases must be linked 
via a continuous (and reversible) path through configuration space. In the context of the method of Mastny and de 
Pablo, this necessitates repeated (back and forth) transitions between the two pure phases. Since no mention is made 
in ref. [17] of any such transitions, the validity of the method would appear to be questionable [19]. 

Another very recent approach to the F-CS problem, which bears some resemblance to that of ref. [17], is the "multi- 
NpH" method of Escobedo [20] . Here a path is constructed in the enthalpy of the system, spanning the range of values 
between those typical of the respective equilibrium phases at some prescribed pressure. The system temperature is 
ascertained at each point along the path, and exhibits a discontinuity at some value of the enthalpy as the favored 
phase of the system changes. TI with respect to the temperatiire variations along the path yields the free energy 
difference. The problem reported with this method was that it was apparently necessary for equilibration purposes to 
initialize the system in the CS phase at all points along the enthalpy path. However, it was difficult to be sure whether, 
for a given enthalpy value, the system had relaxed to the phase of minimum free energy, and hence whether the TI 
result was unbiased. This seems to us to be a manifestation of the familiar drawback of TI, namely the difficulty of 
finding a reversible integration path. 

Phase Switch Monte Carlo (PSMC) [21] is a further relatively recent method designed to directly link F and CS 
phases via a reversible sampling path. The method builds on previous work [22] in which it was demonstrated 
that the configuration spaces associated with two phases of a many-body system can both be visited in a single 
Monte Carlo simulation, by harnessing extended sampling methods to facilitate a direct switch from one phase to the 
other. The method, which is quite generally applicable, was initially deployed in a study of hard sphere freezing [21]. 
Subsequently, it was applied to soft potentials by Errington [23], who used it to calculate two points on the freezing 
line of the Lennard-Jones system for a number of system sizes. 

The present paper complements and extends Errington's PSMC study of the LJ system. It is organized as follows. 
In sec. II we review some of the key results in the statistical mechanical formulation of F and CS phase configurational 
weights, and show how to construct, on this basis, a theoretical framework for computational estimation of free energy 
differences. Sec III describes how this framework is exploited in principle, and realized in practice, by the PSMC 
method. An application of the method to the Lennard-Jones system is reported in sec. IV where we present estimates 
of large portions of the freezing curve for a number of system sizes. Tlic^sc^ we compare with the results of a previous 
GDI study, bootstrapped using TI at a single coexistence point. Finally sec. V details our conclusions and the 
prospects for further applications. 
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II. STATISTICAL MECHANICS 



Within the framework of statistical mechanics, the free energy of a given phase is expressible in terms of its 
configurational weight. Information regarding this weight can be accumulated via a simiilation in the course of an 
exploration (sampling) of the microstates of the phase. While for fluid phases such an exploration encounters no great 
obstacles, a complication arises in the case of crystalline solids. Specifically, it transpires that the relevant phase space 
of the solid is effectively fragmented into a number of mutually inaccessible regions. Each fragment corresponds to a 
distinct permutation of the particles with respect to the lattice sites. The solid phase configurational weight can only 
be estimated for the single fragment in which it is initiated. Since on symmetry grounds the weight contribution of 
every fragment is identical, the overall configurational weight of the solid is obtained by multiplying the measured 
weight of one fragment by the number of fragments. In this section we first show how to correctly count this fragment 
number. We then proceed to describe how the ratio of configurational weights of two phases is related to the ratio of 
their a priori probabilities (a quantity directly accessible to simulation) and thence to the free energy difference. 

Consider a periodic system of N classical particles confined to a volume V, variable under a prescribed constant 
external pressure p, and in equilibrium with a heat bath at prescribed constant temperature T. Within this constant- 
NpT ensemble the configurational weight of a phase may be written as 



Zj{N,p,T)= dVe-P''Z^{N,V,T). (1) 
Jo 

Here, 7 labels the phase, while 

-'V 

Z^{N,V,T) = l[ drie-*(«), (2) 

with $ the configurational energy [24]. The 7-label on the integral stands for some configurational constraint that 
picks out configurations {r} that 'belong' to phase 7. In the present work we shall be concerned with phases which 
are either fluid (7 = F) or crystalline solid (7 = CS) in character, and choose to formulate the constraint in a 
form which reflects the situation actually encountered in Monte Carlo simulations of single phases. Specifically, let 
Rj . . . R^^ = {R}'' denote some representative configuration of phase 7. Then the constraint may be regarded as 
picking out those configurations which can be reached from {R}'^ on the simulation timescale. The timescale is 
presumed to be sufficiently long to explore 'one phase' but still short compared to (unaided) inter-phase traverses. 

It is convenient to use the sites defined by {R}''' as the origins of the particle coordinates. Thus we define a set of 
displacement vectors {u} where 



Ui = ri-R'J , (3) 

and write <I>7({u}) = cl>({R''' + u}). In the case of the fluid phase, particles have the run of the entire system and 
hence all contributing configurations are reachable from any one; we may write simply 



Zf{N, V,T) = r\ rfu,e-*-«"» , (4) 

where {R}^ is some specific but arbitrary fluid configuration, which can be selected at random in the course of MC 
exploration of the fluid phase. 

For the crystalline phase it is natural to choose {R}'^'^ to define the sites of a lattice of the appropriate symmetry 
and scale. But here one must recognize that the complete CS configuration space actually comprises a number 
of distinct mutually inaccessible fragments [25] corresponding essentially to the different permutations of particles 
amongst lattice sites. By symmetry each fragment should contribute equally to the configurational weight. But in 
general, the typical timescale for particle interchanges between lattice sites greatly exceeds the accessible simulation 
timescale. In these circumstances MC simulation will visit (and thus count) only the states within the fragment in 
which it is initiated. The total configurational weight of the CS phase must therefore be obtained by multiplying the 
measured contribution of one fragment by the total number of fragments. Since global translation (permitted by the 
boundary conditions) ensures that one fragment includes all possible locations of any chosen particle, the number of 
fragments is the number of ways of assigning the other — 1 particles to N —1 Wigner-Seitz cells of some underlying 
notional fixed lattice [26]. This number is {N — 1)!. Thus 
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Zcs{N, V,T) = {N-l)\\\ duie-*^-«"» . (5) 

In order to make contact with the simulation methodology to be described below, it is expedient to define the total 
a priori probability of phase 7. This is given by the ratio of the configurational weight of phase 7 to the total weight 
of the two phases [5] : 



nW^,P,-l) zAN,p,T) + Zcs{N,p,T) ■ 

Combining this last equation with eqs. 1, 4 and 5, allows the ratio of the configurational weights of the two phases to 
be expressed in terms of the ratio of their a priori probabilities [22] : 



^ Zf{N,p,T) _ P{F\N,p,T) 
^'^^ Zcs{N,p,T) P{CS\N,p,T) 

{N - 1)! dVe-py nf=i Iv,in}cs du,e-^cs(M) ' 
The link between the configurational weight of a phase and its Gibbs free energy is forged by the definition: 



Z^(iV,p,T) = e-«-(^'f^^). (8) 

Inserting this into eq. 7, the Gibbs free-energy-density difference between the phases can be written in the (strategically 
suggestive) form: 



Ag = gcs{N,p,T) - gF{N,p,T) = ^ ln7eF,cs • (9) 

Eqs. 7 and 9 provide a foundation for MC studies of fluid-solid phase coexistence. Their formal promise is as follows: 
given an inter-phase sampling scheme that visits both the F and CS phases, T^-FjCS (and hence Ag) may be obtained 
as the ratio of the a priori probabilities of the phases. The latter quantity can be estimated simply as the ratio of the 
number of visits to each phase recorded in the course of a simulation. In order to exploit this framework, however, 
one must first design a singularity-free phase space path linking the F and CS phases, and then formulate a sampling 
strategy to traverse it. Phase Switch Monte Carlo is one realization of such a scheme. 



III. PHASE SWITCH MONTE CARLO 



The relative stability of the F and CS phases is determined by the ratio of the associated configurational weights 
^F,CS (Eq. 7). To measure that ratio one needs a MC procedure which visits both solid and fluid regions of con- 
figuration space in the course of a single simulation run. The phase switch MC method is a general approach that 
facilitates such two-phase sampling. Its principal feature is a phase space 'leap' [22] that directly maps a pure phase 
configuration of one phase onto a pure phase configuration of the other. The motivation for choosing such an inter- 
phase route for the F-CS problem is that it circumvents interfacial (mixed-phase) configurations and their attendant 
sampling difficulties (see sec. I). In this section we describe how to apply the phase switch method in this context. 

The key to implementing the phase switch is the representation (eq. 3) of particle coordinates in term of displace- 
ments with respect to a representative configuration {R}'''. By construction, the system may be transformed between 
the CS and F phase representative configurations simply by switching the vectors (Rf ^ Rf'^Vi). Thus by con- 
tinuity, any CS (F) configuration 'sufficiently close' to the representative one will also transform to a F (CS) state 
under this operation [27]. This phase switch can itself be realized as an MC step, so that the phase label 7 becomes 
a stochastic variable. By 'sufficiently close' here, we mean that the energy change associated with the phase switch is 
such as to afford a reasonable chance of acceptance. Fig. 1 provides an illustration of the mechanism. 

We term the set of configurations for which the MC phase switch will be accepted the 'gateway' states. In general, 
however, these gateway states constitute only a small fraction of the total respective configuration spaces of the phases. 
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Fluid Crystalline solid 



FIG. 1: Schematic illustration of the phase switch mechanism. The dots identify the representative sites {R}^ in each phase; the 
displacement vectors {u} connect the centers of the distinguishable (numbered) particles to these sites. The switch operation 
shown swaps the representative sites of one phase for those of the other phase, whilst maintaining {u} constant. More generally, 
any desired transformation of {u} (eg. a scaling -see text) can be incorporated in the switch [5]. The particular configuration 
{u} shown is a "gateway" state (see text) because the magnitude of the effective energy change under the switch is small. 




Fluid Crystalline solid 

FIG. 2: Schematic illustration of the phase switch operation in terms of the regions of configuration space associated with the 
F and CS phases. A bias (dashed arrows) is constructed such as to enhance the probability of the subsets of "gateway" states 
(the white islands) within the single-phase regions, from which the switch operation (the large arrow) will be accepted. Note 
that the switch accesses only one of the (A'^ — 1)! CS phase space replica fragments (see text). 

and consequently will rarely (if ever) be visited on simulation timescales. To ensure effective two-phase sampling, the 
probabilities with which the gateway states are visited must be enhanced. This can be achieved by appeal to extended 
(biased) sampling methods. Fig. 2 shows a schematic representation of the procedure. 

The bias necessary to promote sampling of the gateway states is administered with respect to an 'order parameter' 
-some suitably chosen macrovariable of the system. The freedom in selecting this order parameter is circumscribed 
by the requirement that the associated microstates form a contiguous path through phase space linking the large 
number of equilibrium (typical) microstates to the relatively small number of gateway states. In the next subsection, 
we describe one choice of order parameter which has proved itself successful in this regard. 

A. Order parameter and extended sampling strategy 

One can devise a variety of order parameters that serve to form a suitable path leading from the equilibrium states 
to the gateway states. A definition that we have found to be effective is a variant of one originally employed by one 
of us in ref. [21], and which closely resembles a recent proposal by Errington [23]. Here the order parameter comes 
in two parts (or modes): 'tether' and 'energy'. The tether mode serves to draw particles close to the representative 
sites to which they are nominally associated. Then, once all particles are within a prescribed distance of these sites, 
tether mode switches off and an energy mode order parameter takes over to guide particle to the gateway states for 
which the phase switch energy cost is sufficiently small to be accepted. 

To elaborate, let Mm,-y denote the order parameter in mode m and phase 7. Then for tether mode we write m = T 
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and define an associated order parameter 

1/2 



-^r,7({u}) = I ]y X] ™ax{0, Ui - Uc} j 



(10) 



wliere Ui = juij/y^/"^ is the distance of particle i from its lattice site measured in units of the box length, and Uc is a 
prescribed dimensionless threshold radius. Only particles whose displacement Ui exceeds this threshold contribute to 

The tctlu^r iiiodc' is active iff ui > Uc for at least one particle i, i.e. when Mq-^^ > 0. Otherwise control hands over 
to the 'energy' mode m = £; its associated order parameter is defined by 

M£,^{{u}) = sgn(Afy^)ln(l + \ASy^\) , (11) 

where 

A£y^ = - £f) - {£y - ) (12) 

measures the change (under the phase switch operation) of the Hamiltonian £j{{u.}, V) = $-y({u}) +pV with respect 
to its value £'^^^ in the representative microstate {R}'"', with the latter scaled to the instantaneous volume of phase 7 

[23, 28]. The presence of the logarithm in eq. 11 is designed to moderate the scale of the contribution of the energy 
cost to Ms^-y, which might otherwise become excessively large for particles with a strongly repulsive core to their 
interaction potential. 

Given these definitions, and recalling eqs. 4 and 5, the entire region of accessible configuration space can be described 

by the ensemble 

Z{N,p,T,{n}) = Y, dVj{ duie-«-«">'^) , (13) 
where T-l.y is the effective Hamiltonian defined by 

n^{{u},V) = mu},{R}r)+pV + 7jm,^{M) 

-5^,csin{N-l)\ (14) 

Here {77} represents a set of multicanonical weights [5, 29], associated with the macrostates of the macrovariable 
M — Mfn.'y- The set can usefully be split into four 'branches', one for each combination of mode and phase. The role 
of the weights is to modify (with respect to the Boltzmann sampling) the probability with which the various M values 
are visited. More specifically, as discussed below, we tailor their form in order to ensure roughly uniform sampling 
over the relevant range of M. 

Simulations in the ensemble described by eq. 13 allow one to measure the multicanonical probability distribution 

N 

P{M,V,£,j\N,p,T,{r,}) = Z-'Y[ / duie-«-«">'^)5(M - M({u}, y))5(f^ - f^({u}, V)) , (15) 

i=i •'t 

the form of which is accumulated in practice via a simple list (see sec. Ill B 2). The bias introduced by the multicanon- 
ical weights can be unfolded to give the equilibrium distribution (to within an unspecified normalization constant): 

P{M,V,£,^\N,p,T) = P{M,V,£,'y\N,p,T,{r]})e^'--''^''^ (16) 

The distribution of any single observable can be readily extracted from P{M, V, £, 7) by integration. Similarly the 
ratio of configurational weights is obtained as 
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Zcs{N,p,T) 

J dM dV d£ P{M,V,£,F\N,p,T) 
J dM dV d£ P{M, V, £, CS\N, p, T) 



(17) 



from which the Gibbs free-energy-density difference follows directly via eq. 9. 



B. Implementation 



1. MC moves 



The Monte Carlo procedure wc have adopted utilizes four types of update, which wc describe below. In implementing 
these updates, it is computationally expedient to hold one particle fixed at its representative site in each phase. This 
suppresses the global translation mode in the CS phase and eliminates the need for association swaps (see below) in 
this phase. Choosing to make the assignment Vi^^ = R-7=Af (implying that Ui^^ =0), wc use {u}^, to represent sets 
{u} of displacement coordinates which satisfy this condition. Consequences for acceptance probabilities are examined 
in appendix. A. 

The MC procedure for one-phase exploration has three types of update: 'Particle translations', 'Association swaps' 

and 'Volume moves'. 

1. 'Particle translations'. This is a standard MC procedure: a site (identified by one of the vectors in {R}'-^'^ 
or {R}^) is chosen at random and a candidate state is chosen by incrementing the position coordinate of the 

particle associated with that site by a random vector whose components are drawn from a zero-mean uniform 
distribution of prescribed width. This operation changes both {r} and {u}^, 

2. 'Association swaps'. In this operation we choose two distinct sites at random {i and j say) and identify the 
corresponding displacement vectors Uj and uj. The candidate state is defined by the replacements 



This process can be regarded as a change of association: the particle which was associated with site j is now 
associated with site i (and vice versa). It changes the representation of the configuration (the coordinates {u}); 
but it leaves the physical configuration invariant. 

We apply it only to the fluid phase where particles can wander far from their representative sites and need to 
be reined back in order to reach the gateway states. One may implement association updates in the CS phase 
too: the simulation can then be thought of as exploring different CS- phase fragments; the factor of {N — 1)! in 
Eq. 5 is no longer needed and it is no longer necessary to clamp one particle. 

3. 'Volume moves'. The volume is also updated in the conventional way, by a random walk of prescribed maximum 
step size, with particle position coordinates {r} and representative sites {R}'' rescaled . Note that we maintain 
the ratio Vy /Vj constant throughout, so the notional 'conjugate' phase 7' also undergoes a dilation by the same 
factor. 

In each of the above cases the transition to the candidate state is accepted (see appendix A) with the probability 



Uj — > u • = Uj + Rj - Rj 
Uj — > Uj = Uj -I- Ri - Kj 



(18) 
(19) 



p„({u}.,y^{uK,y' |7) = 

mm{l,exp[-AH^+N\n{V'/V)]} 



(20) 



with 



AH-y = n^{{u}:, V) - n-y{{uu V) . 



(21) 
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All three types of MC move may involve a change in the order parameter M and hence a change A77 in the mul- 
ticanonical weights implicit in Ti (cf. eq. 14). Note that particle translations or association swaps that hand over 
control from one order parameter mode to another, involve a weight change which is calculated from the respective 
branches of the weight function. For example, for a hand-over from tether mode to energy mode, in which the order 
parameter changes from Mr to Ms the appropriate weight change is calculated as 

Ar? = ve,^{M£,^) - r/r.^C^r.^) • (22) 

4 ^Inter-phase switcK . The final type of MC update is the phase switch, which entails replacing one set of the 
representative configuration vectors, {R}''' say, by the other, {R}'' . If the equilibrimn densities of the two 
phases were close to one another it would also be possible to keep the volumes constant in the switch. But that 
is not the case in general. Accordingly -if the procedure is to be efficient- the switch needs to incorporate an 
appropriate volume scaling of the system [30]. It is natural to fix that scaling so that a 'typical' volume of 
phase 7 is matched to a 'typical' volume V-^i of phase 7'. In practice we take to be the mean volume of phase 
7. The switch is accepted (see appendix A) with the probability 

p,(7,{uK,y^7',{u};,y') = 

min|l,exp -A?iy^ + (TV + 1) ln(yy /V:y)j | (23) 

where 



AWy^ = Wy({uK,y') - H-,{{u},,V) (24) 

The phase switch occurs only from states in which the energy mode order parameter is active; the associated 
change in multicanonical weights is 

Ar] = r]£,y{M'£^^,)-r]£,j{M£,^) , (25) 

where ,^,=— M^^^ are the new and old order parameter values respectively. This weight discontinuity is at 
our disposal; we tune it so as to cancel the typical contribution to the acceptance probability associated with 
the switch. 

When discussing our results (sec. IV) we shall occasionally refer to a "Monte Carlo sweep" . This we take to comprise 
N trial particle translations, N trial association swaps, 1 trial volume move and 1 trial phase switch. 

2. Data accumulation and histogram extrapolation 

In the course of the simulations, observations are stored as lists [31] Vj,Mj,£j,^j {j = 1,2,.. .). Estimates for the 
distribution of any observables were accumulated by appropriate binning and weighting of the members of the list 

P(X|7V,p,T)oc5;]e''(^^),5(X-X,), (26) 

j 

with X = V;M: etc. 

In analyzing the data, extensive use was made of histogram extrapolation (HE). This method greatly enhances 
computational efficiency by facilitating the scanning of a region of pressure and temperature around each simulation 

state point, without recourse to further simulation. To achieve this. HE reweights the data obtained at one state 
point {p,T) to yield estimates appropriate to another not-too-distant state point {p',T'). The reliable range of the 
extrapolation in p and T is dependent on the statistical quality of the original data. For further information on the 
HE technique, we refer the interested reader to ref. [33]. Its utility in the context of tracking a freezing line will be 
illustrated in sec. IV. 
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3. Weight function construction 

Prior to collection of phase switch data, it is necessary to determine a suitable set of multicanonical weights {77} 
to enable the system to pass readily between the equilibrium and gateway states. In practice the requisite weights 
are defined with respect to the macrostates formed by discretising (binning) the order parameter macrovariable 
M = Mm,^- For a given mode m and phase 7, the bin width AM was chosen such that for all pairs of successive bins 
M(*) and M(*+i) = M^*) + AM, the corresponding weights satisfy 7]{M^'^+'^^)/r]{M'^''>) < 2. Doing so ensures that 
the weight associated with any given macrostate is representative of the full range of the underlying macrovariable, 
thereby guaranteeing a reasonable MC acceptance rate. 

Consideration of eqs. 14 and eqs. 15 shows that the form of the weights that confers equal probability on every 
macrostate M (and hence ensures that all are well sampled) is r/(A'/) — — In P{M). For adequate sampling one 
therefore requires a prior estimate of P(M). A variety of techniques exist for providing such an estimate, some of 
which are described in refs. [5, 34]. The approach we have utilized in this work is the transition matrix Monte Carlo 
method (TMMC) [23, 34-36]. 

TMMC extracts information on macrostate probabilities by focusing on the statistics of transitions between these 
macrostates. The method proceeds by defining a collection matrix, C, in which information regarding the macrostate 
transitions is accumulated. At every MC step, the proposed transition is recorded in C (irrespective of whether or 
not it is actually accepted) according to: 



C(M^M') ^C(M^M') +a 

C(M^M) ^C(M^M) +(l-a), (27) 

where a is the acceptance probability of the move under the Hamiltonian £j{{u},V) (rather than the effective 
Hamiltonian W^). 

The macrostate transition probabilities derive from the collection matrix via: 



T(M^M') = ^ , (28) 

EfcC(M^Mfe)' ^ ' 

where the denominator on the RHS sums over all possible values of the macrostates. Knowledge of these transition 
probabilities permits, in turn, an estimate of the probabilities P{M). Although a number of approaches exist (see 
eg. ref. [23]), for extracting P{M) from T, a simple, yet unbiased and stable method, considers exclusively those 
transitions which occur between adjacent macrostates (i.e. examines only the diagonal and first off-diagonal elements 
of T). Application of the 'balance' condition for macrostate transitions [34]: 



P{M') ^ TjM^M') 

P{M) T(M'^M) ' ^ 
then permits assignment of the weight function via 

,(M')-,(M) = -ln(2^ga^,). (30) 

The weight function was updated in this manner at regular intervals (typically every 1000 Monte Carlo sweeps). 

Each update allowed the simulation to explore an ever wider range of order parameter values [23]. We note, however, 
that weight function updates should not be performed too frequently during the weight evolution process because, 
strictly speaking, they violates detailed balance. Empirically, however, at the stated update interval, the solution was 
found to be highly accurate and the procedure efficient. 

For the purposes of constructing the weight function, we found it beneficial to decompose the full range of M into 
a number of overlapping 'windows'. Within each window a fragment of the overall weight function was determined 
and joined for continuity to those of neighboring windows. The weight differences across the three mode transitions 
(tether to energy mode for each phase, as well as the phase switch itself) were found using a simple root-finding 
algorithm in a simulation restricted to the appropriate extremes of the order parameter. Both energy modes were 
calculated as single windows, while the solid and fiuid tether modes were split into 3 and 10 windows respectively (for 
the largest system size studied N = 500). 
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IV. APPLICATION TO THE LENNARD-JONES SYSTEM 



A. Model and strategy 



We have employed the Phase Switch Monte Carlo method to study the freezing properties of particles whose 
interactions are described by the Lennard-Jones potential 



The system is assumed to be contained within a cubic box of volume V, which fluctuates under the control of an applied 
external pressure p and is in thermal equilibrium at temperature T. In common with previous studies of freezing in 
the LJ system [11, 23], particle interactions were truncated at one-half the box length, and running tail corrections 
applied to the energy in the standard manner [3]. As described in sec. IIIB 1, the MC moves employed were particle 
translations, volume moves, association swaps and phase switches. Of these, we note that volume moves and phase 
switches are computationally inexpensive to implement because the associated energy change is obtainable directly 
from the scaling with (respect to coordinate dilations) of the repulsive and attractive parts of the configurational 
enc;rgy [37]. Association swaps are similarly cheap because they leave the configurational energy unaltered. Since the 
position of one particle was clamped at its representative site in each phase, CS phase association swaps were not 
required (see sec. IIIBl and appendix A). 

The phase switch MC method requires for its operation the specification of a representative configuration {R}''' 
for both the CS and F phases. For {R}'^^ we adopt a perfect lattice of the appropriate symmetry (here fee) and 
scale (corresponding to a rough estimate of the CS phase coexistence density). For {R}^, we randomly select a 
configuration in the course of a simulation of the fluid phase [32] . 

With regard to the order parameter (sec sec III A), the hand-over from tether to energy mode is controlled by the 
parameter Uc in eq. 10. This should be chosen such that -at the hand-over- particles are sufficiently close to their 
representative sites that the energy cost A£~^/~^ of the switch (cq. 12) is not saturated at a large value. Otherwise Mg^-^ 
will be unable to guide the particles to positions favorable for the phase switch. We find that satisfactory results are 
obtained by choosing Uc consistent with Mg^^ « 5 at the hand-over. 

The formalism described in Sees. II and III assumes that the configurational constraint confines the CS phase to 
the phase space fragment in which it is initiated. While this appears to be strictly observed for F-CS coexistence in 
hard spheres, as studied in ref. [21], for the softer interactions of the LJ system, interchange of neighboring particles 
between lattice sites was observed to occur, albeit very rarely. The simplest way of dealing with this effect is to 
suppress it; to which end we introduced an upper bound on the particle displacements in the CS phase. We set 
jwij'"'^ < 0.65(7, a choice which (we have verified) has a negligible eSect on the equilibrium properties of the CS phase 
for the system sizes we have studied. 

The strategy for obtaining the F-CS coexistence line is as follows. Initially one requires a rough estimate for some 
coexistence state point. This may derive either from any available literature estimates or, alternatively, coexistence 
can be bracketed by scanning (backwards and forwards) a range of pressures at some nominated temperature until the 
system spontaneously transforms between phases. Next it is necessary to customize a suitable set of multicanonical 
weights which allow the system to sample both phases. This can be done by employing, for example, the transition 
matrix method described in sec. IIIB 3. An example weight function is shown in fig. 3 for a system oi N = 256 
particles close to coexistence at /3 = e/kT = 0.6. 

Having determined a suitable weight function, a production run is then performed to obtain data of high statistical 
quality. From this data one can readily extract (see sec. IIIB 2) the ratio 7?,f,cs (eq- 17) and hence the free energy 
difference Ag (eq. 9) between the phases. Histogram extrapolation (HE) (sec. IIIB 2) then permits the location of the 
coexistence pressure at which this difference vanishes (i.e. when 7^f,CS = !)• The corresponding coexistence number 
densities can then be simply read off from the peak positions in the coexistence form of P{p), where p = N/V. Fig. 4 
shows a typical coexistence form of P{p), obtained for N = 256 at /3 = 0.6. 

In addition to simplifying the determination of the coexistence pressure at a given temperature, HE also permits 
simultaneous extrapolation in the temperature. This facilitates the scanning of a segment of the coexistence line 
in the p — T plane. Although the range of near coexistence p — T values for which a single simulation provides 
reliable information is necessarily limited, additional simulations can be performed near the vestiges of this range. 
The data from simulations performed at a number of different state points can then be synthesised self-consistently 
using multiple histogram reweighting [33] . As well as guiding the choice of near coexistence values of p and T for 
subsequent simulations, HE also provides an estimate of the appropriate order parameter distribution P{Mm,-f), which 
(recall sec. Ill B 3) serves as a suitable weight function. Thus by a repeated combination of simulation followed by 
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Bin number n 

FIG. 3: An example weight function used for a near-coexistence production run. The data shown arc for a system of A'^ = 256 
LJ particles at /3 = 0.6, P = 13.7. The four branches of the weight function have been laid out with respect to the underlying 
discretisation into bins (see sec. Ill B 3) . 



25- 




FIG. 4: The coexistence form of the density distribution P{p) for a system of iV = 256 LJ particles at /9 = 0.6, p = 13.722. A 
selection of data points are shown. Lines are guides to the eye. 



HE, the coexistence line can be straightforwardly tracked, without the need to ever re-determine a weight function 
ah initio [31]. 

B. Results 

Estimates of the freezing boundary have been obtained using the PSMC method within the constant-iVpT ensemble 
for systems of sizes N = 108, 256 and 500 particles. Freezing boundary data has also been obtained for N = 32, but 
here the system is sufficiently small that transitions back and forth between F and CS phases occur spontaneously, 
over a range of pressures, and a density distribution (sampling both phases) can be determined and a coexistence 
pressure inferred- without the need for PSMC. 

Spontaneous transitions were also found to occur (albeit much more rarely) in the N = 108 system. Here we 
observed spontaneous freezing as we tracked the coexistence curve using the methods described in sec. IV A. We 
emphasize that this effect was also seen for this system size in a 'bare' constant- A^pT simulation at coexistence and 
does not therefore appear to be caused (or rendered more frequent) by the multicanonical weighting of PSMC. The 
larger system sizes {N = 256 and N = 500) did not exhibit spontaneous transitions on the timescale of our simulations, 
presumably because of the increased free energy cost of forming an interface. 

The spontaneous freezing that occurs in the A^ = 108 system undermines our strategy whereby the phase label 7 
keeps track of the current phase. Nevertheless, we were able to suppress the problem by tracking the coexistence 
line along a path parallel to its true trajectory, but displaced slightly towards the fluid side. This was achieved by 
simulating at a pressure some 5% less than the coexistence value predicted by the histogram extrapolation from the 
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previous ncar-cocxistcncc state point. The estimate of the coexistence pressure was easily corrected for the imposed 
shift by means of histogram extrapolation. 
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FIG. 5: (a) The estimated F-CS phase diagram in the pressure-inverse temperature plane for the four systems sizes studied 
in this work. The data shown derive from 20 separate simulation state points for the N = 32 system size, 37 state points 
for A'^ = 108, 17 points for N = 256 and 4 points for the A'^ = 500 system size. Also included for comparison are the PSMC 
estimates of Errington [23] for /? = 4/3, 1/2, together with the GDI results of Agrawal and Kofke [11] for = 500. (b) A 
closeup of the region around = 4/3. The vertical error bars corresponding to Errington's data points [23] and the horizontal 
error bars to the GDI study. Symbols are the estimates arising from the present study (given explicitly in tables I-III), with 
uncertainties smaller than the symbol size in each case. Lines are interpolations between the data points, based on histogram 
extrapolation. 

The measured phase boundary in the p — T plane, and a portion of the phase diagram in the p — T plane, are 
shown in figs. 5 and 6 respectively for the various system sizes studied. Full data are also tabulated in tables I-III 
together with their uncertainties. With regard to the latter, we note that analysis of statistical errors within the 
PSMC framework is extremely transparent: at a given temperature, the uncertainty in the pressure is simply given 

by 



= ' ^^^^ 

where /S.V = Vf — Vcs, and cr[7^F,Cs] is the associated uncertainty in 7?.f,cS: this being controlled (at heart) by the 
statistics of the inter-phase switch and measurable from a simple block error analysis. We note that owing to the 
extensive factor AV in the denominator of eq. 32, the minimum uncertainty in T^-p^cs required to attain a prescribed 
error in the pressure grows like N [38] . 
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FIG. 6: A portion of the phase diagram in the p — T plane in the region of the triple point. Shown are the estimated solid 
and fluid coexistence densities for A'^ = 500 and N = 108. Also included for comparison are the GDI estimates of Agrawal and 
Kofke[ll] for N = 500. Uncertainties are smaller than the symbol sizes; lines are guides to the eye 

As with any method, the error analysis becomes more involved if one chooses to combine a number of separate data 

sets via a multiple histogram analysis [33]. In these circumstances we have found it useful to perform a bootstrap 
error analysis with 100 re-samples, taking an estimate of the block size as the correlation length. The error estimate 
was found to depend only weakly on the block size. All error bars in tables I-III are calculated in this way and 
correspond to a 67% confidence interval. We stress that it is not strictly necessary to perform this more sophisticated 
error analysis, adequate (albeit slightly overly conservative) error bars can be assigned simply on the basis of block 
averages. 

Ideally one should like to extrapolate coexistence parameters obtained for a range of finite system sizes to the 
thermodynamic limit using a finite-size scaling ansatz. In the hard sphere case, a clear dependence of the 
coexistence pressure has been observed using the PSMC method [21. 23], and there seems no reason to expect a 
different scaling in the case of the LJ potential. However, in the present study (as well as in that of ref. [23]), a good 
fit to this scaling form could not be obtained. As noted in ref. [23], the reason for this is presumably traceable to 
the truncation of the interparticle potential cutoff, which was set to be one-half the box length in all cases. This 
choice, while facilitating direct comparison with existing literature data (principally ref. [11]), necessarily introduces 
a degree of coupling between the potential truncation and the system size, notwithstanding the use of a mean-field 
based tail correction. For the rather limited range of system sizes studied in this work, the use of such a tail correction 
is unlikely to provide a good approximation to the full potential, particularly for the smaller system sizes. Indeed, it 
is well known that many aspects of the phase behaviour of the LJ system are acutely sensitive to the details of the 
potential cutoff [23, 28, 40]. 

V. DISCUSSION AND CONCLUSIONS 

In summary, the Phase Switch MC method allows one to locate directly and transparently fluid-solid coexistence 
parameters and their associated uncertainties within the appropriate (constant-pressure) ensemble. To achieve this 
the method connects the phases via a direct inter-phase sampling path, thereby facilitating estimates of free energy 
differences from a single simulation. The course of this path encompasses configurations which are exclusively pure- 
phase in character. Accordingly, one can be confident that it robust, i.e. free of both singularities and crystalline 
defects. 

We have applied the PSMC method to calculate, for a number of system sizes, significant portions of the p — T 
freezing curve of the LJ fluid to high accuracy (see the error bars quoted in tables I-III). In so doing we have 
illustrated the strategic advantages of histogram extrapolation (HE) as an efficient means of tracking the coexistence 
line. Furthermore we have seen (cf. refs. [13, 31]) that HE provides reliable estimates of the requisite multicanonical 
weight function, thereby obviating the need to redetermine the weights from scratch at each successive near-coexistence 
state point. While it is probably fair to say that the PSMC remains a somewhat computationally intensive strategy 
(at least in the context of the F-CS problem), it is by no means prohibitively so on the scale of its competitors. Given 
the strengths of the method we have identified, we feel it should represent an attractive option, especially when high 
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TABLE I: Solid-fluid coexistence curve data for the = 108 system. Tabulated in columns 1 — 6, respectively, are the inverse 
freezing temperature, pressure, solid number density, fluid number density, energy per particle (solid) and energy per particle 
(fluid). Numbers in parentheses indicate the the 67% confidence limit for the rightmost digit(s). Note that due to the steepness 
of the coexistence curve at low pressures (cf. fig. 5), it is expedient in this regime to determine the coexistence temperature at 
prescribed pressure. 

precision estimates of freezing boundaries are required. 

Wc have compared our results with those of the previous PSMC study of Errington [23] and the GDI study of 
Agrawal and Kofke [11] for TV = 500. Not surprisingly, our results agree (to within error) with those of Errington. 
However, (and as also noted by Errington), the results for A'' = 500 do not overlap with the error bars quoted in the 
GDI study of the same system size (and potential truncation) by Agrawal and Kofke [11]. This is possibly due to the 
inherent limitations of the GDI scheme as already noted in sec. I -specifically the lack of any means to reconnect the 
phases beyond the initial 'starting' coexistence state point from where the; integration is launched. It may additionally 
reflect a poor estimate of the starting point itself, which was determined using TI by transforming the L J system to 
a hard sphere reference system via a route which takes in an intermediate system of soft spheres. In particular, the 
coexistence pressure taken for the hard sphere reference system is now believed to be an overestimate [21, 23, 41]. 
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TABLE II: Solid-fluid coexistence curve data for A'^ = 256. See caption of tab. I for details. 



AT = 500 

e/kT pa^/e pes (^^ pF Scs/Ne Ep/Ne 

1.336(2) 0.5500 0.9655(1) 0.8566(2) -6.637(1) -5.486(2) 

1.3613(3) 0.3836 0.96337(4) 0.85274(3) -6.815(1) -5.6705(2) 

1.393(1) 0.2091 0.9617(1) 0.8490(1) -7.007(1) -5.871(1) 

1.418(2) 0.0200 0.9581(1) 0.843(1) -7.198(1) -6.08(1) 

TABLE III: Solid-fluid coexistence curve data for N = 500. See caption of tab. I for details. 

As regards the prospects for further applications, it would certainly be worthwhile in the specific case of the LJ 
system to attempt to decouple the system size scaling of the cutoff from intrinsic finite-size effects by repeating the 
present study for a constant interparticlc potential cutoff. This should permit a reliable extrapolation of coexistence 
parameters to the thermodynamic limit, albeit for the particular choice of cutoff. More generally, it would be 
interesting to attempt to apply the PSMC method to molecular systems. Here, however, the order parameter would 
probably need to be augmented in order to not only draw molecules to their representative sites, but also to align 
them appropriately. 
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APPENDIX A: ACCEPTANCE CRITERIA 

Here we derive the Monte Carlo acceptance probabilities used in the PSMC method. 

As described in the sec. HI, we choose to fix one particle at its representative site. The effect of this is to reduce 
the configurational weight of each phase by a factor of the system volume V. Thus for the fluid phase, eq. 4 becomes 
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while for the CS phase (using eq. 5), we have 
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In analogy to eq. 13 the phase space relevant to the problem is described by the multicanonical ensemble 
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The configurational probability follows as 
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Consider now a MC move (particle translation, association swaps, or volume move) which leaves the phase label 
unchanged. Let 
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be the probability that the system is found in the region of configuration space of volume Y\- duidV around {u}*, V. 
Similarly let 



JV-l 



(A7) 



represent the probability associated with the range of configuration space reached by implementing a change 
({u}*,y) — > ({u}^,y). Detailed balance requires 



PajA^B) ^ Pb 

Pa{B ^ A) Pa 



(A8) 



where 



Pb _ V eM-n,i{u}i,V'mti dWjdy 
Pa Vexp[-H^{{u},,V)]Y\fj;^du,dV 

Now for a volume update that increments the current value, dV — dV, while 



V 



N 



(A9) 



(AlO) 
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Thus 

PajA^B) _ exp[~n^{{u}i,V')+N\nV'] 
Pa{B ^ A) exp[-W^({u}^, V) +N\nV] 

from which the acceptance probability follows as 



p„({u}.,y^{u};,y' |7) = 

mm{l,exp[-AHj+N\n{V'/V)]} (All) 

with 

AH^ = H^{{u}l V) - W^({u}„ V) . (A12) 

Wc note that the acceptance formula eq. All is actually no different to that pertaining to a fully unconstrained 
constant- A^pT ensemble [3]. Indeed a little thought reveals that the act of fixing one particle is simply equivalent to 
viewing the unconstrained system from the reference frame of that particle. Such a change of reference frame has no 
consequences for acceptance criteria. 

Consider next the inter-phase switch. This replaces one set of representative vectors, {R}''' say, by the other, {R}'^ , 
while at the same time scaling the system volume and particle displacements by a constant factor chosen such that a 
'typical' volume Vj of phase 7 is matched to a 'typical' volume Vj' of phase 7'. For a volume update which scales the 
current value, one has dV'/dV = V'/V, and referring back to the detailed balance relation eq. A9 one finds 



'Yl. 
V 



N+l 



(A13) 



The switch acceptance probability follows as 



p„(7,{uK,V^7',{uK,n = 

min{l,exp AH^j + {N + l)ln{VY /V-y) } 



(A14) 



where 



AHjy = ({u}:,, V) - n-y{{uU V) . (A15) 

Note that a different formulation of the switch acceptance probability applies when the displacement vectors are held 
constant during the switch operation [21, 23, 30]. 
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